Bienvenid@s a la primera tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Videos de las clases:
Documentación:
Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.
Respuesta
Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.
Respuesta
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:
# Manipulación de estructuras
library(tidyverse)
library(dplyr)
library(tidyr)
# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
# Manipulación de varios plots en una imagen.
library(gridExtra)
# Análisis bayesiano
library("rethinking")
Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:
install.packages("rethinking")
En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:
install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')
El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.
En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:
\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \] donde \(\Gamma(\alpha)\) es una constante, usualmente se le llama función gamma.
De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).
Respuesta
# Algoritmo Metropolis-Hastings
metropolis_hastings <- function(theta0, alpha, beta, iterations) {
gamma_density <- function(x, alpha, beta) {
if (x > 0) {
return(x^(alpha - 1) * exp(-beta * x))
} else {
return(0)
}
}
samples <- numeric(iterations)
samples[1] <- theta0
for (t in 2:iterations) {
theta_star <- rnorm(1, mean = samples[t - 1], sd = 1)
if (theta_star > 0) {
acceptance_ratio <- gamma_density(theta_star, alpha, beta) /
gamma_density(samples[t - 1], alpha, beta)
} else {
acceptance_ratio <- 0
}
u <- runif(1)
if (u < acceptance_ratio) {
samples[t] <- theta_star
} else {
samples[t] <- samples[t - 1]
}
}
return(samples)
}
# Configuración inicial
alpha <- 5
beta <- 15
theta0 <- 1
# Diferentes cantidades de repeticiones
samples_1k <- metropolis_hastings(theta0, alpha, beta, 1000)
samples_10k <- metropolis_hastings(theta0, alpha, beta, 10000)
samples_100k <- metropolis_hastings(theta0, alpha, beta, 100000)
# Histogramas
par(mfrow = c(1, 3))
hist(samples_1k, breaks = 30, probability = TRUE,
main = "N = 1000", xlab = "Theta", col = "skyblue")
hist(samples_10k, breaks = 30, probability = TRUE,
main = "N = 10000", xlab = "Theta", col = "skyblue")
hist(samples_100k, breaks = 30, probability = TRUE,
main = "N = 100000", xlab = "Theta", col = "skyblue")
A medida que aumenta la cantidad de repeticiones, la distribución obtenida se vuelve más estable y representativa de la distribución objetivo. Con pocas iteraciones (N=1000, N=1000), la muestra puede ser ruidosa e imprecisa, con una mayor variabilidad debido a la dependencia inicial de las propuestas y aceptaciones. Sin embargo, al incrementar el número de iteraciones (N= 100000, N=100000), la muestra converge hacia la distribución gamma real, capturando su forma y reduciendo la influencia de las condiciones iniciales.
# Efecto de la condición inicial
samples_theta0_0.1 <- metropolis_hastings(0.1, alpha, beta, 100000)
samples_theta0_10 <- metropolis_hastings(10, alpha, beta, 100000)
# Histogramas
par(mfrow = c(1, 3))
hist(samples_100k, breaks = 30, probability = TRUE,
main = "Theta0 = 1", xlab = "Theta", col = "skyblue")
hist(samples_theta0_0.1, breaks = 30, probability = TRUE,
main = "Theta0 = 0.1", xlab = "Theta", col = "skyblue")
hist(samples_theta0_10, breaks = 30, probability = TRUE,
main = "Theta0 = 10", xlab = "Theta", col = "skyblue")
La condición inicial (\(\theta_0\)) tiene un impacto significativo en las primeras iteraciones del algoritmo, especialmente si \(\theta_0\) está lejos de los valores más probables de la distribución objetivo. Sin embargo, con un número suficientemente grande de repeticiones, el efecto de \(\theta_0\) desaparece gracias a la convergencia de la cadena de Markov. Esto se evidencia en los experimentos realizados: para \(\theta_0 = 0.1\), \(\theta_0 = 1\), y \(\theta_0 = 10\), todas las cadenas convergen a la misma distribución final cuando se ejecutan \(N = 100000\) iteraciones. Por lo tanto, la condición inicial es importante a corto plazo, pero no afecta el resultado final si el algoritmo se ejecuta por suficientes iteraciones.
# Comparación con la distribución gamma real
samples_real <- rgamma(100000, shape = alpha, rate = beta)
# Graficar comparación
par(mfrow = c(1, 1))
hist(samples_100k, breaks = 30, probability = TRUE,
main = "Comparación con Gamma real", xlab = "Theta", col = "skyblue")
lines(density(samples_real), col = "red", lwd = 2)
legend("topright", legend = c("MH Samples", "Real Gamma"),
col = c("skyblue", "red"), lwd = 2)
Se observa que el histograma obtenido con el algoritmo Metropolis-Hastings para \(N = 100000\) coincide de manera consistente con la densidad de la distribución gamma. Esto valida que la implementación del algoritmo es correcta y que las muestras generadas se distribuyen según la distribución objetivo con un número suficiente de iteraciones.
A work by CC6104